RTCH  6 


REPORT  DOCUMENTATION  PAGE 


form  ADDrovea 
QMS  Nc  0704^01  SF 


1  AGENCV  USE  ONLY  (Leave  DiariK:  2.  REPORT  DATE 

8/16/96 


4.  TITLE  AND  SUBTITLE 

I  Monte  Carlo  Aialysis  of  Nonequilibrium  Ionized  Flows 


I  6.  AUTHOR(S) 

Iain  D.  Boyd 


3.  REPORT  TYPE  AND  DATES  COVERED 

Final  Technical  Report  7/1/94  to  6/30/96 


5.  FUNDING  NUMBERS  ! 

,  PE  -  61102F 

.d  Flows  _  2308  ! 

SA  -  AS  I 

-  G  -  F49620-94-1-0328  i 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 
Cornell  University 
246  Upson  Hall 
•  Ithaca,  NY  14853 


AFOSR-TR-96 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADORESS(ES) 
AFOSR/NA  4 

110  Duncan  Avenue,  Suite  B115  A  /  JL 

Bolling  AFB,  DC  20332-0001  ///  K 


10.  SPONSORING  MONITORING  i 

AGENCY  REPORT  NUMBER  , 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approved  for  public  release; 
distribution  is  unlimited 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

This  final  report  summarizes  the  research  carried  out  by  our  research  group  concerned 
with  the  development  of  a  numerical  approach  for  simulating  the  flows  through  and 
expanding  from  low-power  hydrogen  arcjets.  In  Section  I,  we  describe  the  development 
of  a  particle-based  Monte  Carlo  approach  for  simulating  the  necessaryj^ollision  and 
plasma  effects  that  occur  in  these  electric  propulsion  devices.  In  Section  H,  we  describe 
an  extensive  study  to  make  comparison  between  our  model  predictions  and  existing 
experimental  data.  It  is  concluded  that  the  Monte  Carlo  approach  provides  an  accurate 
means  for  simulating  the  flow  inside  of  the  nozzle  of  the  arcjet  offering  the  potential  to 
improve  arcjet  performance.  It  is  also  concluded  that  accurate  Monte  Carlo  simulation  of 
the  arcjet  plumes  offers  the  potential  to  assess  accurately  integration  of  arcjets  with  real 
spacecraft  geometries. 


14.  SUBJECT  TERMS 


Arcjet  (electric)  propulsion;  Numerical  Simulation 


IS.  NUMBER  OF  PAGES 

43 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  118.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 

OF  REPORT  I  OF  THIS  PAGE  OF  ABSTRACT 


S0-5;00 


S*t‘i'"Od''c  293  ''Rev  j-39- 


TfflS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DTTC  CONTAINED 
A  SIGNIFICANT  NUMBER  OF 
PAGES  WHICH  DO  NOT 
REPRODUCE  LEGffiLY. 


4* 


FINAL  TECHNICAL  REPORT 

For  research  supported  by 
AFOSR  Contract  No.  F49620-94- 1-0328  (Year  2) 

for  the  period  07/01/94  to  06/30/96 

MONTE  CARLO  ANALYSIS  OF  NONEQUILIBRIUM  IONIZED  FLOWS 

prepared  by 
Iain  D.  Boyd  (l) 

Xiaoming  Liu  (2)  and  Jitendra  Balakrishnan  (2) 


School  of  Mechanical  and  Aerospace  Engineering 
Cornell  University, 

246  Upson  Hall 
Ithaca,  NY  14853. 


Work  supported  by 

Air  Force  Office  of  Scientific  Research 
Program  Manager:  Dr.  Mitat  A.  Birkan 

August  1996 


(f)  Principal  Investigator 
(2)  Graduate  Research  Assistant 


19961015  013 


1 


Table  of  Contents 


Page 


Abstract  3 

Overview  4 

1.  Monte  Carlo  Model  of  Hydrogen  Arcjet  Flows  5 

n.  Extensive  Validation  of  the  Monte  Carlo  Model  24 


2 


Abstract 


This  final  report  summarizes  the  research  carried  out  by  our  research  group  concerned 
with  the  development  of  a  numerical  approach  for  simulating  the  flows  through  and  expanding 
from  low-power  hydrogen  arcjets.  In  Section  I,  we  describe  the  development  of  a  particle-based 
Monte  Carlo  approach  for  simulating  the  necessary  collision  and  plasma  effects  that  occur  in 
these  electric  propulsion  devices.  In  Section  II,  we  describe  an  extensive  study  to  make 
comparison  between  our  model  predictions  and  existing  experimental  data,  it  is  concluded  that 
the  Monte  Carlo  approach  provides  an  accurate  means  for  simulating  the  flow  inside  of  the 
nozzle  of  the  arcjet  offering  the  potential  to  improve  arcjet  performance.  It  is  also  concluded  that 
accurate  Monte  Carlo  simulation  of  the  arcjet  plumes  offers  the  potential  to  assess  accurately 
integration  of  arcjets  with  real  spacecraft  geometries. 


Overview 


Arcjets  are  being  installed  on  spacecraft  for  a  variety  of  orbit  maintenance  tasks.  There  is 
a  requirement  for  development  of  physically  accurate  numerical  models  to  describe  the  nozzle 
and  plume  flows  of  these  thrusters.  Accurate  simulation  of  nozzle  flows  offers  the  potential  to 
use  the  model  to  optimize  the  flow  conditions  and  geometry  of  the  arcjet  to  improve  propulsion 
efficiency  which  is  usually  of  the  order  of  30-40%.  This  is  particularly  important  as  the  need  for 
small-scale  thrusters  becomes  an  important  component  in  the  development  of  micro-  and  nano¬ 
satellites.  Accurate  simulation  of  arcjet  plumes  is  necessary  for  prediction  of  spacecraft 
interaction  effects  such  as  contamination  of  solar  arrays,  torques,  and  surface  heating. 

The  modeling  effort  reported  here  has  led  to  the  development  of  a  Monte  Carlo  model  for 
describing  the  nozzle  and  plume  flows  of  a  1  kW  arcjet.  This  work  focused  on  the  flow  of 
hydrogen  through  the  arcjet  due  to  the  availability  of  a  large  volume  of  detailed  experimental 
data  for  this  configuration.  Development  of  the  Monte  Carlo  model  required  implementation  of 
high-temperature  kinetics  and  plasma  models  for  the  hydrogen  system.  This  work  is  described  in 
Section  I.  The  Monte  Carlo  approach  offers  the  possibility  of  simulating  nonequilibrium 
transport  and  chemical  kinetics  effects  compared  to  continuum  formulations.  In  addition,  direct 
simulation  of  the  important  plasma  effect  of  ohmic  heating  is  included  in  the  Monte  Carlo 
technique  for  the  first  time.  Ihe  inclusion  of  ohmic  heating  is  required  to  predict  the  specific 
impulse  measured  experimentally.  This  work  was  presented  as  an  invited  paper  (AIAA  96-2022) 
at  Ae  AIAA  27th  Fluid  Dynamics  Conference,  New  Orleans,  Louisiana,  June  1996. 

An  extensive  series  of  comparisons  between  the  Monte  Carlo  predictions  and  detailed 
experimental  data  was  conducted  in:  (1)  the  interior  nozzle  flow;  (2)  the  nozzle  exit  plane;  and 
(3)  the  arcjet  plume.  The  comparisons  of  properties  in  the  arcjet  plume  are  the  first  of  their  kind. 
Continuum  numerical  methods  are  generally  unable  to  compute  the  plume  flows  due  to  their  very 
low  density.  These  comparisons  are  described  in  Section  II.  Generally,  very  good  agreement 
between  computations  and  measurements  was  obtained  for  flow  properties  such  as  number 
densities,  velocities,  and  temperatures  for  the  interior  nozzle  flow,  for  the  nozzle  exit  plane,  and 
for  the  plume  expansion.  Good  agreement  was  also  obtained  with  measurements  of  specific 
impulse.  The  only  flow  field  property  for  which  there  was  significant  disagreement  was  the 
number  density  of  atomic  hydrogen  which  was  consistently  predicted  to  be  a  factor  of  3  to  4 
lower  than  values  measured  experimentally  at  the  Phillips  Laboratory.  This  work  was  presented 
as  a  contributed  paper  (AIAA  96-3297)  at  the  32nd  AIAA/ASME/SAE/ASEE  Joint  Propulsion 
Conference,  Orlando,  Florida,  July  1996.  An  additional  paper  presented  at  this  conference 
(Wysong,  Pobst  and  Boyd,  "Comparisons  of  Hydrogen  Atom  Measurements  in  an  Arcjet  Plume 
With  DSMC  Predictions,"  AIAA  Paper  96-3185)  also  resulted  from  this  research  effort. 

The  conclusion  drawn  from  these  studies  is  that  the  Monte  Carlo  approach  is  an 
appropriate  method  for  accurate  modeling  of  the  nozzle  and  plume  flows  of  hydrogen  arcjets. 
The  arcjets  manufactured  by  Olin  Aerospace  that  are  being  flown  in  space  use  hydrazine  as 
propellant.  To  model  the  flows  of  the  real  thrusters,  extension  of  the  Monte  Carlo  model  is 
required  to  include  the  more  complicated  kinetics  mechanisms  of  the  H2-N2  system. 
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I.  MONTE  CARLO  MODEL  OF  HYDROGEN  ARCJET  FLOWS 


Abstract 

A  particle-based  Monte  Carlo  numerical  method  is  developed  for  computation  of  flow  in  a  low-thrust 
hydrogen  arcjet.  The  method  employs  the  direct  simulation  Monte  Carlo  technique  to  compute  the  nonequi¬ 
librium  fluid  mechanics  and  thermochemical  relaxation.  Simulation  of  plasma  effects  is  included  at  a  simple 
level  that  ensures  charge  neutrality.  A  new  model  is  developed  to  include  ohmic  heating  in  the  simulation. 
Results  are  presented  for  a  flow  condition  that  has  been  studied  experimentally  with  a  number  of  different 
diagnostic  techniques.  A  strong  degree  of  thermal  nonequilibrium  is  found  in  the  flow.  It  is  demonstrated 
that  ohmic  heating  is  important,  and  that  the  Monte  Carlo  model  captures  the  important  physics  very  well. 
It  is  also  found  that  the  simple  model  employed  for  motion  of  the  charged  species  is  inadequate  leading  to 
poor  prediction  of  some  flow  properties. 


Introduction 

Arcjets  operated  at  a  power  level  of  about  1  kW  are  an  attractive  candidate  for  use  as  electric  propulsion 
devices  on  spacecraft.  They  offer  substantial  increases  in  specific  impulse  compared  with  chemical  rockets 
and  resistojets.  Unfortunately,  the  performance  of  arcjets  remains  poor  with  propulsion  efficiencies  in  the 
range  of  30  to  50%.  There  is  a  need  for  accurate  numerical  simulation  of  the  flows  generated  by  these  devices 
both  for  improving  their  performance,  and  for  analysis  of  potential  spacecraft  interaction  effects  such  as 
contamination.  There  are  many  coupled  physical  phenomena  at  work  in  an  arcjet  that  makes  numerical 
simulation  diflacult.  Nonequilibrium  kinetic  processes  include  rotational  and  vibrational  relaxation,  and 
several  chemical  mechanisms  including  dissociation,  recombination,  and  ionization.  In  addition  plasma 
effects  are  important. 

A  hydrazine  arcjet  has  been  developed  and  flown  in  space.^  In  terms  of  basic  research,  however,  hydrogen 
arcjets  have  received  considerable  attention.  Several  experimental  studies  have  been  conducted  on  an  arc-jet 
designed  and  fabricated  at  NASA  Lewis  Research  Center.  A  schematic  diagram  of  this  thruster  is  shown 
in  Fig.  1.  The  NASA  Lewis  hydrogen  arcjet  has  been  used  in  experimental  investigations  to  determine 
its  performance^  and  flow  flelds.^“^  Of  considerable  importance  is  the  fact  that  several  different  diagnostic 
techniques  have  been  applied  to  these  flows  yielding  measurements  that  include  translational  and  rotational 
temperature,  velocity,  and  number  densities  of  atomic  and  molecular  hydrogen.  These  data  offer  excellent 
opportunities  for  detailed  calibration  of  modeling  efforts.  Arcjet  models  based  on  continuum  formulations 
have  been  developed  by  several  groups.^"®  There  are  two  main  problems  with  the  continuum  approach. 
First,  there  is  a  significant  degree  of  thermal  and  chemical  nonequilibrium  in  these  flows.  Second,  it  is  very 
difficult  to  extend  the  continuum  simulations  out  to  the  plume  produced  by  the  arcjet.  This  means  that 
these  methods  cannot  be  used  to  predict  spacecraft  interference  effects. 

This  paper  describes  the  development  of  the  direct  simulation  Monte  Carlo  method  (DSMC)  for  appli¬ 
cation  to  a  hydrogen  arcjet.  In  satellite  propulsion,  the  DSMC  technique  has  been  applied  and  successfully 
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verified  to  the  plume  of  a  hydrazine  rocket®,  and  to  the  nozzle  and  plume  flows  of  both  a  nitrogen  resisto- 
jet^®  and  a  cold  flow  of  hydrogen  through  an  arcjet.^^  Preliminary  DSMC  results  for  heated  arcjet  flows  were 
reported  in  Ref.  5.  In  this  paper,  description  is  provided  of  the  collision  and  chemistry  models  developed  for 
this  application.  In  addition,  the  implementation  of  plasma  effects  are  discussed.  In  particular,  this  study 
describes  a  new  model  for  computation  of  ohmic  heating  using  the  DSMC  technique.  Results  are  presented 
to  illustrate  the  general  nonequilibrium  nature  of  the  arcjet  flows  and  comparison  is  made  with  experimental 
data. 


Numerical  Method 

The  direct  simulation  Monte  Carlo  method  (DSMC)  uses  the  motions  and  collisions  of  particles  to 
perform  a  direct  simulation  of  nonequilibrium  gas  dynamics.  Each  particle  has  coordinates  in  physical 
space,  three  velocity  components,  and  internal  energies.  A  computational  grid  is  employed  to  group  together 
particles  that  are  likely  to  collide.  Collision  selection  is  based  on  a  probability  model  developed  from  basic 
concepts  in  kinetic  theory.  The  method  is  widely  used  for  rarefied  nonequilibrium  conditions  and  finds 
application  in  hypersonics,  materials  processing,  and  micro-machine  flows. 

The  extension  of  the  technique  for  application  to  arcjet  flows  requires  development  of  a  number  of  new 
physical  models.  Arcjets  are  characterized  by  high  temperatures  (as  high  as  30,000  K)  which  give  rise  to 
strongly  ionized  flows  (as  much  as  30%  ionized).  In  terms  of  intermolecular  collisions,  mechanisms  must 
be  included  in  the  numerical  model  to  account  for:  (1)  momentum  transfer;  (2)  rotational  excitation;  (3) 
vibrational  excitation;  (4)  dissociation/recombination;  and  (5)  ionization/recombination.  The  relatively  high 
level  of  ionization  requires  modeling  of  plasma  effects  at  some  level.  Each  of  these  modeling  issues  is  discussed 
in  the  following. 

Momentum  Transfer 

For  a  hydrogen  arcjet,  the  four  main  species  present  are  H2,  H,  and  e"^.  In  general,  the  most 
important  collision  types  that  determine  the  arcjet  flow  field  are  H2-H2,  H2-H,  H-H,  H2“e*”,  and  H-e"*.  In 
the  present  work,  the  momentum  transfer  cross-sections  for  the  first  three  collision  types  are  obtained  from 
the  work  of  Vanderslice  et  al.^^  in  which  collision  integrals  were  tabulated  for  the  temperature  range  of  1,000 
to  15,000  K.  These  collision  integrals  are  first  converted  into  total  collision  cross  sections  and  are  shown  in 
Fig.  2.  Using  this  data,  appropriate  parameters  are  determined  for  the  Variable  Hard  Sphere  (VHS)  collision 
model^^  that  is  widely  used  in  DSMC  computations.  In  the  present  DSMC  code,  a  single  value  of  the  viscosity 
temperature  exponent  cj  that  appears  in  the  VHS  model  is  used  for  all  collision  types.  This  places  a  small 
restriction  on  the  level  of  agreement  that  is  achievable  for  all  collisions.  The  VHS  parameters  identified 
as  providing  the  best  fit  to  the  data  are:  a;=0.82,  T7.g/=300  K,  with  reference  collision  cross-sections  of 
(7^2-^2=3.554  X  10"^®  m^,  o-i/2_7/=3.140  x  10"^®  m^,  and  (7//_i/=3.554  x  10"^®  m^.  Comparison  of  the 
VHS  models  to  the  data  of  Ref.  12  is  made  in  Fig.  2.  The  viscosities  obtained  using  these  VHS  parameters 
are  compared  with  the  computations  of  Vanderslice  et  al.  as  a  function  of  temperature  in  Fig.  3.  The  level  of 
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agreement  is  quite  good.  Also  shown  in  Fig.  3  is  the  viscosity  obtained  when  VHS  parameters  for  hydrogen 
for  low  temperature  flows  are  employed.  In  this  case,  the  VHS  parameters  are:  a;=0.67,  Tre/=273  K,  <th2-H2 
=  o'H2-h  =  o-h-h  =  2.606  X  10”“^®  m^.  These  parameters  were  used  in  the  previous  study  that  applied  the 
DSMC  technique  to  cold  hydrogen  flow  through  an  arcjet  geometry.^^  These  results  provide  an  indication  of 
the  sensitivity  of  the  flow  field  modeling  to  momentum  transfer. 

The  VHS  reference  cross-section  for  the  H2-e'“  collision  is  determined  to  be  cr//2-e“  =3.547  x  10“^^  m^ 
using  the  experimental  data  of  Crompton  and  Sutton. Comparison  of  the  model  and  the  data  is  shown 
in  Fig.  4.  The  range  of  relative  velocity  computed  for  the  VHS  model  corresponds  to  the  thermal  speed  of 
electrons  over  the  temperature  range  from  1,000  to  15,000  K.  Clearly,  the  VHS  model  does  not  capture  the 
data  very  well.  This  represents  a  requirement  for  further  development  in  the  future.  The  VHS  reference 
cross-section  for  the  H-e“  collision  is  determined  to  be  =7.118  x  10”^^  m^  using  the  simple  analytical 

expression  given  by  Yos.^^  Comparison  of  the  VHS  model  with  this  expression  is  shown  in  Fig.  5  and  reveals 
generally  very  good  agreement  except  at  very  low  temperatures. 

Other  cross-sections  are  much  less  important  and  are  assigned  values  such  that  the  results  of  hydrogen 
atoms  are  used  for  hydrogen  ions.  For  example,  the  cross-sections  shown  in  Fig.  2  for  H2-H  collisions  are 
also  used  for  H2-H'^  collisions.  While  this  results  in  some  lack  of  accuracy,  it  is  expected  that  the  effect  will 
be  small  as  these  collisions  are  relatively  infrequent. 

Rotational  and  Vibrational  Relaxation 

Molecular  hydrogen  is  an  interesting  diatomic  molecule  in  that  it  has  relatively  large  values  for  the 
characteristic  temperatures  of  rotation  and  vibration.  This  results  in  long  relaxation  times  for  these  modes. 
In  the  present  work,  the  transfer  of  energy  between  the  translational  and  rotational  energy  modes  employs 
the  phenomenological  approach  of  Borgnakke-Larsen.^®  In  this  scheme,  a  probability  of  exchanging  energy 
is  applied  to  each  collision.  If  it  is  decided  that  energy  transfer  does  occur,  the  post-collision  energies 
are  sampled  from  equilibrium  energy  distributions.  Previously,  a  probability  model  for  determining  the 
probability  of  translational-rotational  energy  transfer  for  hydrogen  at  low  temperatures  was  developed  and 
verified. This  same  model  is  employed  again  here  although  it  is  recognized  that  it  may  be  inaccurate  at 
the  very  high  temperatures  encountered  in  the  arcjet. 

Vibrational  energy  is  modeled  using  the  quantized  oscillator  form  of  the  Borgnakke-Larsen  approach 
developed  by  Boyd  and  Bergemann.^^  The  probability  of  vibrational  energy  exchange  employs  the  famil¬ 
iar  Landau-Teller  relaxation  time  with  parameters  obtained  experimentally  by  Kiefer.^®  Once  again  it  is 
recognized  that  this  approach  may  be  inaccurate  at  the  very  high  temperatures  found  in  the  arcjet. 

The  average  probabilities  of  rotational  (RT)  and  vibrational  (VT)  energy  exchange  for  H2-H2  collisions 
are  compared  as  a  function  of  temperature  in  Fig.  6. 

Chemical  Reactions 

There  are  only  two  chemical  reactions  implemented  in  the  present  study.  They  are  dissociation  of 
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molecular  hydrogen: 


and  ionization  of  atomic  hydrogen: 


H^M  i — ^iI++e“+M 


where  in  each  case  M  is  any  available  species.  The  dissociation  rate  for  M=H2,  H,  and  is  obtained  from 
Ref.  19  as  1.82  x  10""^^  rj^-0,251  -51,275/T)  m^/sec.  For  dissociation  through  electron  impact,  the  rate 

is  taken  from  McCay  and  Dexter^^  as  3.17  x  10~“^®  exp(  -51,275/T)  m^/sec.  The  ionization  rate  for 
M-H2,  H,  and  H+  is  taken  from  Ref.  20  as  7.00  x  10"^^  T"*^  °  exp(  -157,825/T)  m^/sec.  For  electron  impact 
ionization,  Ref.  20  gives  the  rate  as  5.27  x  10^  T"^’°  exp(  -157,825/T)  m^/sec.  It  is  not  possible  to  capture 
this  temperature  dependence  using  the  standard  DSMC  chemistry  model.  In  particular,  the  temperature 
exponent  for  this  collision  pair  must  be  no  smaller  than  about  -1.  An  expression  that  provides  reasonable 
agreement  with  the  actual  rate  is  1.32  x  10"^  T"^  °  exp(  -120,000/T)  m^/sec.  Comparison  of  the  exact  rate 
and  the  expression  used  in  the  DSMC  computation  is  shown  in  Fig.  7. 

The  dissociation  and  ionization  reactions  are  modeled  using  the  DSMC  chemistry  model  developed  by 
Boyd.^^  In  the  case  of  dissociation,  the  coupling  to  the  vibrational  energy  is  set  to  zero.  Recombination 
reaction  rates  are  computed  from  the  above  forward  rates  and  the  equilibrium  constants.  These  reactions 
are  then  simulated  using  the  DSMC  recombination  model  described  in  Ref.  21. 


Plasma  Modeling 

In  principle,  it  is  possible  to  model  plasma  effects  using  the  particles  of  the  DSMC  algorithm  through 
the  addition  of  Particle  In  Cell  (PIC)^^  steps  into  the  computer  code.  In  practise,  there  is  such  a  large 
discrepancy  between  the  plasma  and  the  collision  length  scales  that  need  to  be  resolved  in  the  arcjet  flow, 
that  this  approach  is  infeasible.  Instead,  a  very  simple  scheme  is  employed  in  which  the  electrons  are 
constrained  to  move  through  the  flow  field  with  the  average  ion  velocity  in  each  cell.  This  ensures  charge 
neutrality,  but  obviously  omits  the  charge  field  effects. 

A  second  important  plasma  effect  that  is  included  in  the  DSMC  simulation  is  ohmic  heating.  This  is  an 
important  source  of  energy  that  results  from  the  current  represented  by  the  flux  of  charged  species  through 
the  thruster.  The  electrical  resistivity  is  given  by: 


7]  = 


me  Et  ^ei 
riee^ 


(1) 


where  e  is  the  charge  of  an  electron,  me  is  the  electron  mass,  rig  is  the  electron  number  density,  and  Ud  is  the 
mean  collision  rate  between  an  electron  and  heavy  particle  i.  Power  generated  in  the  plasma  due  to  ohmic 
heating  is: 

n  =  7?j2.  (2) 
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Here,  J  is  the  current  vector  so  that 


—  e^Y^niUi  -  rieUef'  +  [UiVi  -  rieVe)^] 


(3) 


where  Wj  is  the  ion  number  density,  Ui  and  Vi  are  the  mean  ion  velocities  in  the  z  and  r  directions,  and  Ue 
and  Ve  are  the  mean  electron  velocities  in  the  z  and  r  directions.  The  mean  collision  rates  required  by  Eq.  (1) 
for  e“-H2  and  e“-H  interactions  are  computed  directly  in  the  DSMC  simulation.  For  the  purposes  of  ohmic 
heating,  a  more  accurate  value  of  the  e“-H+  interaction  is  computed  using  the  Coulombic  cross-section: 


<  <^e--H+  >- 


jlogiA) 


247r(eofcTe)2 

where  Co  is  the  permittivity  of  free  space,  Te  is  the  electron  temperature,  and 

127r{€okTe)^ 


A  = 


Q  — 

e^Ue  2 


(4) 


(5) 


The  total  ohmic  power  is  calculated  in  each  cell  of  the  simulation  using  Eqs.  (1-5).  The  energy  repre¬ 
sented  by  n  is  added  to  the  thermal  energies  of  the  charged  species  (electrons  and  ions)  in  each  cell.  The 
average  energy  to  be  added  to  each  charged  particle  in  the  cell  per  time-step  is  then  given  by: 

nAt 


E  = 


nAV 


(6) 


where  At  is  the  simulation  time  step,  AV  is  the  cell  volume,  and  n  is  the  total  number  density  of  charged 
particles.  The  energy  addition  is  performed  such  that  momentum  is  conserved  on  average.  This  is  achieved 
by  increasing  only  the  thermal  velocity  components  of  each  charged  particle  by  an  amount  proportional  to 

Ve. 


Simulation  Details 

The  collision  and  plasma  models  are  implemented  in  an  existing  DSMC  code  that  is  designed  for  efficient 
performance  on  vector  supercomputers.^^  Modest  changes  have  also  been  made  to  allow  the  code  to  execute 
efficiently  on  work-stations.  The  addition  of  the  ohmic  heating  model  adds  only  a  20%  overhead  to  the 
performance  of  the  DSMC  code. 

The  DSMC  computations  are  begun  at  a  location  downstream  of  the  nozzle  throat.  The  initial  profiles 
of  flow  properties  are  obtained  from  continuum  simulations  provided  by  the  code  of  Butler.^  The  interaction 
of  the  gas  with  the  wall  of  the  nozzle  is  assumed  to  be  diffuse  with  full  energy  accommodation  to  the  wall 
temperature  profile  obtained  from  the  continuum  simulation.  Charged  species  are  neutralized  upon  contact 
with  the  wall.  The  computations  extend  through  the  nozzle  exit  plane  to  a  distance  of  30  mm  from  the 
thruster  along  the  axis.  Expansion  is  assumed  to  occur  into  a  perfect  vacuum. 

The  grids  employed  typically  consist  of  500  by  55  cells.  Due  to  the  wide  dynamic  range  in  density  for 
these  flows,  the  cells  are  adapted  to  the  local  mean  free  path,  and  the  time  step  in  each  cell  is  scaled  by  the 
local  mean  time  between  collisions.  Despite  the  time  step  variation,  the  transient  phase  of  the  simulation 
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remains  long.  Typically,  it  requires  at  least  15,000  iterations  before  steady  state  is  reached.  Macroscopic 
results  are  then  obtained  by  sampling  properties  over  a  further  10,000  time  steps.  At  steady  state,  there  are 
at  least  300,000  particles  present  in  the  simulation. 

Results 

In  terms  of  results,  this  investigation  focuses  on  some  of  the  nonequilibrium  aspects  of  hydrogen  arcjet 
flows.  More  detailed  comparisons  of  experimental  data  with  results  obtained  using  the  DSMC  code  described 
in  this  paper  are  reported  in  Refs.  24-25.  The  case  considered  here  is  for  a  power  level  of  P=1.4  kW,  and  a 
mass  flow  rate,  ?7i=13.3  mg/sec.  These  are  conditions  that  have  been  investigated  experimentally  in  a  series 
of  investigations  by  Cappelli  and  his  co-workers  at  Stanford  University,  and  at  the  Phillips  Laboratory. 
First,  results  are  presented  that  are  obtained  using  the  cross-sections  of  Vanderslice  et  al..  Next,  the  effect 
of  the  ohmic  heating  model  is  assessed.  Finally,  sensitivity  of  the  solutions  to  changing  the  cross-sections  for 
neutral  interactions  is  considered. 

Standard  Computation 

The  following  results  are  obtained  using  the  collision  and  plasma  models  described  earlier.  In  Fig.  8a, 
the  concentrations  of  molecular  and  atomic  hydrogen  across  the  nozzle  exit  plane  are  shown.  Experimental 
measurement  of  H2  was  performed  by  Beattie  and  Cappelli^®  using  Raman  scattering.  The  measurements  for 
H  are  those  reported  by  Wysong  et  al..^^  In  general,  the  agreement  between  computation  and  measurement 
is  better  for  H2  than  for  H.  In  particular,  the  simulations  predict  an  atomic  concentration  on  the  axis  that 
is  a  factor  of  3  to  4  smaller  than  the  experimental  data. 

In  Fig.  8b,  the  translational  temperature  profile  of  atomic  hydrogen  across  the  nozzle  exit  plane  is 
shown.  Two  sets  of  experimental  data  reported  by  Wysong  et  al.^^  and  by  Storm  and  Cappelli^^  are  shown. 
In  general,  the  DSMC  results  show  a  small  underprediction  of  the  temperature,  although  the  agreement  is 
quite  good.  The  degree  of  thermal  nonequilibrium  of  the  flow  in  the  nozzle  exit  is  illustrated  in  Fig.  8c. 
Here,  the  translational,  rotational,  and  vibrational  temperatures  of  molecular  hydrogen  predicted  by  the 
DSMC  technique  are  compared  with  the  rotational  temperature  measurements  of  Beattie  and  Cappelli.^® 
The  predicted  rotational  temperatures  offer  satisfactory  agreement  with  the  measured  data.  The  peak 
vibrational  temperature  is  about  a  factor  of  two  higher  than  the  peak  translational  temperature. 

A  further  indication  of  the  degree  of  thermal  nonequilibrium  in  this  flow  is  illustrated  in  Fig.  8d  in 
which  the  electron  temperature  in  the  nozzle  exit  plane  is  shown.  The  experimental  data  is  due  to  Storm 
and  Cappelli.^®  The  DSMC  results  show  a  relatively  high  degree  of  statistical  scatter  due  to  the  small  number 
of  particles  in  the  simulation  that  represent  the  electrons.  Nevertheless,  it  is  clear  that  the  behavior  close 
to  the  axis  is  very  well  captured  by  the  simulation.  Note  that  the  peak  electron  temperature  is  an  order  of 
magnitude  higher  than  the  peak  translational  temperature.  In  Fig.  8e,  the  electron  concentration  across  the 
nozzle  exit  plane  is  shown.  The  experimental  data  was  also  reported  in  Ref.  28.  The  agreement  between 
simulation  and  measurement  is  very  poor  with  the  DSMC  technique  predicting  concentrations  that  are  an 
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order  of  magnitude  higher  than  the  data. 

Finally,  the  axial  component  of  velocity  in  the  nozzle  exit  plane  is  shown  in  Fig.  8f.  The  experimental 
data  is  that  of  Wysong  et  al..^^  The  DSMC  results  are  consistently  lower  than  the  measured  data  across 
the  entire  exit  plane.  In  Fig.  8g,  the  axial  velocity  along  the  nozzle  centerline  is  shown.  The  experimental 
data  is  that  of  Storm  and  Cappelli.^®  Note  that  these  data  were  taken  at  a  slightly  higher  arcjet  power  of 
P=1.48  kW.  It  is  evident  that  the  DSMC  simulation  predicts  lower  velocities  than  the  measured  data  along 
the  entire  extent  of  the  nozzle.  It  should  also  be  noted  that  the  specific  impulse  predicted  by  the  DSMC  code 
is  780  sec  in  comparison  to  a  value  of  about  900  measured  at  NASA  Lewis  Research  Center  under  similar 
conditions.^ 

Effect  of  Ohmic  Heating 

A  second  simulation  is  performed  in  which  the  ohmic  heating  model  is  not  employed  in  the  DSMC 
computation.  It  is  found  that  the  absence  of  ohmic  heating  most  affects  the  temperatures  in  the  computed 
flow  field.  In  Fig.  9a,  the  translational  temperature  of  atomic  hydrogen  in  the  nozzle  exit  plane  is  shown.  By 
comparison  with  Fig.  8b,  it  is  clear  that  the  peak  translational  temperature  is  reduced  by  about  600  K  when 
there  is  no  ohmic  heating.  Clearly,  these  DSMC  results  offer  much  poorer  agreement  with  the  data  compared 
to  those  obtained  with  the  ohmic  heating  model.  In  Fig.  9b,  the  translational,  rotational,  and  vibrational 
temperatures  predicted  in  this  case  are  compared  with  the  measurements  for  rotational  temperature.  The 
vibrational  temperature  is  slightly  higher  in  this  case  due  to  slower  relaxation  caused  by  the  reduction  in 
translational  temperature.  The  rotational  temperature  is  lower  than  that  obtained  with  the  ohmic  heating 
model.  Once  again,  the  agreement  with  the  data  obtained  without  ohmic  heating  appears  to  be  worse  than 
that  shown  in  Fig.  8c.  The  profile  of  electron  temperature  in  the  nozzle  exit  plane  is  compared  with  the 
measured  data  in  Fig.  9c.  In  this  case,  the  peak  electron  temperature  predicted  by  the  simulation  is  only 
about  2,000  K  and  is  clearly  much  lower  than  the  values  observed  experimentally. 

Effect  of  Viscosity 

A  final  simulation  is  performed  in  which  the  momentum  cross-sections  for  interaction  between  the  neutral 
atoms  and  molecules  are  those  employed  in  the  earlier  DSMC  study  of  a  cold  flow  hydrogen  arcjet. The 
ohmic  heating  model  is  again  employed.  In  general,  the  DSMC  results  are  found  to  be  relatively  insensitive 
to  viscosity.  The  most  pronounced  difference  arises  in  velocity.  In  Fig.  10,  the  variation  of  axial  velocity 
along  the  nozzle  for  these  cross-sections  is  compared  with  the  experimental  data.  With  this  viscosity  model, 
there  are  more  collisions  simulated  and  hence  the  gas  accelerates  to  a  slightly  higher  velocity  offering  better 
correspondence  to  the  data.  It  is  believed  that  this  improved  agreement  is  fortuitous,  although  it  is  recognized 
that  the  momentum  cross-sections  of  Ref.  12  may  not  be  accurate  throughout  the  arcjet  as  they  are  only 
evaluated  to  a  temperature  of  15,000  K  whereas  the  peak  temperature  in  the  arcjet  flow  is  about  30,000  K. 

Discussion  and  Conclusions 
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The  first  conclusion  to  be  drawn  from  the  comparisons  of  DSMC  predictions  and  experimental  measure¬ 
ments  is  that  the  model  for  ohmic  heating  appears  to  represent  the  plasma  physics  very  well.  The  increases 
in  the  translational,  rotational,  and  electron  temperatures  obtained  with  the  ohmic  heating  model  all  match 
very  well  to  the  experimental  data.  Also,  in  terms  of  the  thermal  state  of  the  gas,  the  DSMC  results  indicate 
that  a  strong  degree  of  thermal  nonequilibrium  exists  in  the  nozzle  exit  plane  of  the  thruster.  It  is  therefore 
a  recommendation  that  continuum  methods  developed  for  application  to  arcjet  nozzle  flows  should  include 
finite  rates  of  relaxation  for  the  vibrational  energy  of  molecular  hydrogen,  and  a  separate  energy  equation 
for  the  translational  motion  of  the  electrons.  These  are  important  loss  mechanisms  for  thruster  performance 
and  may  explain  why  continuum  methods  that  omit  these  processes  tend  to  predict  higher  values  for  specific 
impulse  in  comparison  to  experimental  data.^ 

It  is  also  clear  from  the  present  investigation  that  the  simple  approximation  that  ensures  charge  neu¬ 
trality  is  insuflRcient.  It  is  hypothesized  that  this  approximation  is  responsible  for  the  poor  prediction  by  the 
numerical  method  of  both  velocity  and  electron  concentration.  Inside  the  arcjet,  the  electrons  will  have  very 
large  thermal  speeds  that  would  separate  them  rapidly  from  the  ions.  However,  electrostatic  forces  tend  to 
pull  the  electrons  back  towards  the  ions,  and  at  the  same  time  accelerate  the  ions  towards  the  electrons. 
This  physical  phenomena  is  clearly  not  captured  in  sufficient  detail  by  the  present  approach  of  moving  the 
electrons  with  the  average  ion  velocity.  This  completely  omits  the  particle  acceleration.  Simulation  of  this 
physical  behavior  would  result  in  the  electrons  moving  through  the  flow  at  significantly  higher  velocities  that 
would  in  turn  decrease  the  electron  number  density  throughout  the  nozzle  flow.  Collision  of  the  accelerated 
ions  and  electrons  with  neutral  particles  would  also  result  in  small  increases  in  the  fluid  velocity  throughout 
the  nozzle  flow.  In  addition,  it  is  anticipated  that  such  an  approach  would  also  lead  to  an  increase  in  the 
specific  impulse  computed  by  the  simulation  method. 

As  stated  in  the  earlier  discussions,  the  main  difficulty  associated  with  direct  inclusion  of  the  electric 
field  effects  is  the  large  disparity  in  length  and  time  scales  between  plasma  and  collision  processes.  In  the 
present  investigation,  the  computational  grid  is  sized  according  to  the  collision  length  scales.  For  the  electron 
number  densities  found  at  the  nozzle  throat  of  the  arcjet,  the  cells  there  would  have  to  be  decreased  by  a 
factor  of  1000  to  adequately  resolve  the  plasma  behavior.  This  represents  a  problem  therefore  that  cannot 
simply  be  resolved  using  supercomputers. 
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Fig.  1.  Schematic  diagram  of  the  1  kW  hydrogen  arcjet. 
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Fig.  2.  Average  momentum  transfer  cross-sections  as  a  function  of  temperature. 
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Fig.  3.  Viscosity  coefficient  as  a  function  of  temperature. 


Fig.  4.  Momentum  transfer  cross-sections  for  H2"e  collisions  as  a  function  of  relative  velocity. 
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Fig.  5.  Momentum  transfer  cross-sections  for  H-e  collisions  as  a  function  of  temperature. 


Fig.  6.  Average  probabilities  of  RT  and  VT  energy  exchange  for  H2-’H2  collisions  as  a  function 
of  temperature. 
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Fig.  8a.  Concentrations  of  neutral  species  in  the  nozzle  exit  plane. 


Fig.  8b.  Translation  temperature  of  atomic  hydrogen  in  the  nozzle  exit  plane. 
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Fig.  8g.  Axial  velocity  component  along  the  nozi 
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Fig.  9a.  Translation  temperature  of  atomic  hydrogen  in  the  nozzle  exit  plane:  no  ohmic 
heating. 


Fig.  9b.  Thermal  nonequilibrium  of  molecular  hydrogen  in  the  nozzle  exit  plane:  no  ohmic 
heating. 
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Fig.  9c.  Electron  temperature  in  the  nozzle  exit  plane:  no  ohmic  heating. 


Fig.  10.  Axial  velocity  component  along  the  nozzle  centerline:  neutral  viscosity  model  of 
Ref.  11. 
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II.  EXTENSIVE  VALIDATION  OF  THE  MONTE  CARLO  MODEL 


Abstract 

A  particle-based  Monte  Carlo  numerical  method  is  applied  to  model  the  nozzle  and  plume  flows  of  a 
low-thrust  hydrogen  arcjet.  The  method  employs  the  direct  simulation  Monte  Carlo  technique  to  compute 
the  nonequilibrium  fluid  mechanics  and  thermochemical  relaxation.  Simulation  of  plasma  field  effects  is 
included  at  a  simple  level  that  ensures  charge  neutrality.  A  model  is  implemented  to  include  ohmic  heating 
in  the  simulation.  Results  are  presented  for  two  different  arcjet  flow  conditions.  Comparisons  of  Monte  Carlo 
results  with  experimental  measurements  of  a  variety  of  flow  field  properties  are  made  in  the  nozzle  interior,  in 
the  nozzle  exit  plane,  and  in  the  plume  flow  field.  Comparison  of  prediction  with  measurement  is  also  made 
for  thruster  performance.  In  general,  good  agreement  is  obtained  between  simulation  and  measurement. 
This  indicates  that  the  Monte  Carlo  model  is  a  viable  approach  for  optimization  of  arcjet  performance,  and 
for  estimation  of  spacecraft  interaction  effects. 


Introduction 

Arcjets  operated  at  a  power  level  of  about  1  kW  are  an  attractive  candidate  for  use  as  electric  propulsion 
devices  on  spacecraft.  They  offer  substantial  increases  in  specific  impulse  compared  with  chemical  rockets 
and  resistojets.  Unfortunately,  the  performance  of  arcjets  remains  poor  with  propulsion  efficiencies  in  the 
range  of  30  to  50%.  There  is  a  need  for  accurate  numerical  simulation  of  the  flows  generated  by  these  devices 
both  for  improving  their  performance,  and  for  analysis  of  potential  spacecraft  interaction  effects  such  as 
contamination.  There  are  many  coupled  physical  phenomena  at  work  in  an  arcjet  that  makes  numerical 
simulation  difficult.  Nonequilibrium  kinetic  processes  include  rotational  and  vibrational  relaxation,  and 
several  chemical  mechanisms  including  dissociation,  recombination,  and  ionization.  In  addition  plasma 
effects  are  important. 

A  hydrazine  arcjet  has  been  developed  and  flown  in  space.^  In  terms  of  basic  research,  however,  hydrogen 
arcjets  have  received  considerable  attention.  Several  experimental  studies  have  been  conducted  on  an  arc-jet 
designed  and  fabricated  at  NASA  Lewis  Research  Center.  A  schematic  diagram  of  this  thruster  is  shown 
in  Fig.  1.  The  NASA  Lewis  hydrogen  arcjet  has  been  used  in  experimental  investigations  to  determine 
its  performance^  and  flow  fields.^"^  Of  considerable  importance  is  the  fact  that  several  different  diagnostic 
techniques  have  been  applied  to  these  flows  yielding  measurements  that  include  translational  and  rotational 
temperature,  velocity,  and  number  densities  of  atomic  and  molecular  hydrogen.  These  data  offer  excellent 
opportunities  for  detailed  calibration  of  modeling  efforts.  Arcjet  models  based  on  continuum  formulations 
have  been  developed  by  several  groups.^"®  There  are  two  main  problems  with  the  continuum  approach. 
First,  there  is  a  significant  degree  of  thermal  and  chemical  nonequilibrium  in  these  flows  that  may  not  be 
simulated  accurately  using  a  continuum  formulation.  Second,  it  is  very  difficult  to  extend  the  continuum 
simulations  out  to  the  plume  produced  by  the  arcjet.  This  means  that  these  methods  cannot  be  used  to 
predict  spacecraft  interference  effects. 
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This  paper  describes  the  application  of  the  direct  simulation  Monte  Carlo  method  (DSMC)  for  compu¬ 
tation  of  the  flow  fields  of  a  hydrogen  arcjet.  In  satellite  propulsion,  the  DSMC  technique  has  been  applied 
and  successfully  verified  to  the  plume  of  a  hydrazine  rocket^,  and  to  the  nozzle  and  plume  flows  of  both 
a  nitrogen  resisto-jet^®  and  a  cold  flow  of  hydrogen  through  an  arcjet. In  this  paper,  brief  description  is 
provided  of  the  kinetic  and  plasma  models  developed  for  this  application.  More  details  may  be  found  in 
Ref.  12.  Extensive  comparison  of  the  Monte  Carlo  calculations  with  available  experimental  data  is  made  in 
the  nozzle  interior,  in  the  nozzle  exit  plane,  and  in  the  plume  expansion  flow  fields.  Comparison  is  also  made 
with  measurement  for  thruster  performance. 

Numerical  Method 

The  direct  simulation  Monte  Carlo  method  (DSMC)  uses  the  motions  and  collisions  of  particles  to 
perform  a  direct  simulation  of  nonequilibrium  gas  dynamics.  Each  particle  has  coordinates  in  physical 
space,  three  velocity  components,  and  internal  energies,  A  computational  grid  is  employed  to  group  together 
particles  that  are  likely  to  collide.  Collision  selection  is  based  on  a  probability  model  developed  from  basic 
concepts  in  kinetic  theory.  The  method  is  widely  used  for  rarefied  nonequilibrium  conditions  and  finds 
application  in  hypersonics,  materials  processing,  and  micro-machine  flows. 

The  extension  of  the  technique  for  application  to  arcjet  flows  is  described  in  detail  by  Boyd.^^  This 
required  development  of  a  number  of  new  physical  models.  Arcjets  are  characterized  by  high  temperatures 
(as  high  as  30,000  K)  that  give  rise  to  strongly  ionized  flows  (as  much  £ts  30%  ionized).  In  terms  of  inter- 
molecular  collisions,  mechanisms  must  be  included  in  the  numerical  model  to  account  for:  (1)  momentum 
transfer;  (2)  rotational  excitation;  (3)  vibrational  excitation;  (4)  dissociation/recombination;  and  (5)  ioniza¬ 
tion/recombination.  The  relatively  large  degree  of  ionization  requires  modeling  of  plasma  effects  at  some 
level. 

For  a  hydrogen  arcjet,  the  four  main  species  present  are  H2,  H,  H"^,  and  e”.  In  general,  the  most 
important  collision  types  that  determine  the  arcjet  flow  field  are  H2-H2,  H2-H,  H-H,  H2“e“,  and  H-e".  In 
the  present  work,  the  momentum  transfer  cross-sections  for  the  first  three  collision  types  are  obtained  from 
the  work  of  Vanderslice  et  al.^^  in  which  collision  integrals  were  tabulated  for  the  temperature  range  of 
1,000  to  15,000  K.  The  cross-sections  for  the  H2“e"  collision  is  determined  using  the  experimental  data  of 
Crompton  and  Sutton. The  cross-section  for  the  H-e“  collision  is  determined  using  the  simple  analytical 
expression  given  by  Yos.^^  The  level  of  agreement  between  the  models  used  in  the  DSMC  code  and  the 
original  sources  is  very  good  with  the  exception  of  the  H2-e“‘  collision.  Other  cross-sections  are  much  less 
important  and  are  assigned  values  such  that  the  results  of  hydrogen  atoms  are  used  for  hydrogen  ions.  For 
example,  the  cross-sections  for  H2-H  collisions  are  also  used  for  H2-H+  collisions.  While  this  results  in  some 
lack  of  accuracy,  it  is  expected  that  the  effect  will  be  small  as  these  collisions  are  relatively  infrequent. 

Molecular  hydrogen  is  an  interesting  diatomic  molecule  in  that  it  has  relatively  large  values  for  the 
characteristic  temperatures  of  rotation  and  vibration.  This  results  in  long  relaxation  times  for  these  modes. 
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Previously,  a  rate  model  for  determining  the  probability  of  translational-rotational  energy  transfer  for  hy¬ 
drogen  at  low  temperatures  was  developed  and  verified. This  same  model  is  employed  again  here  although 
it  is  recognized  that  it  may  be  inaccurate  at  the  very  high  temperatures  encountered  in  the  arcjet.  The 
probability  of  vibrational  energy  exchange  employs  the  familiar  Landau-Teller  relaxation  time  with  param¬ 
eters  obtained  experimentally  by  Kiefer  and  Lutz.^®  Once  again  it  is  recognized  that  this  approach  may 
be  inaccurate  at  the  very  high  temperatures  found  in  the  arcjet.  There  are  only  two  chemical  reactions 
implemented  in  the  present  study.  They  are  dissociation  of  molecular  hydrogen  and  ionization  of  atomic 
hydrogen.  The  rate  coefficients  are  obtained  from  Refs,  17  and  18. 

A  very  simple  scheme  is  employed  to  simulate  electric  field  effects.  The  electrons  are  constrained  to 
move  through  the  flow  field  with  the  average  ion  velocity  in  each  cell.  This  ensures  charge  neutrality,  but 
obviously  omits  the  electric  field  effects.  A  second  important  plasma  effect  that  is  included  in  the  DSMC 
simulation  is  ohmic  heating.  This  is  an  important  source  of  energy  that  results  from  the  current  represented 
by  the  flux  of  charged  species  through  the  thruster.  Using  standard  plasma  physics  theory,  a  source  of 
energy  is  computed  in  each  cell  of  the  simulation  due  to  ohmic  heating.  This  energy  is  added  to  the  charged 
particles  in  each  cell  in  each  iteration  of  the  calculation.  The  scheme  is  described  in  detail  in  Ref.  12. 

The  collision  and  plasma  models  are  implemented  in  an  existing  DSMC  code  that  is  designed  for  efficient 
performance  on  vector  supercomputers.^®  The  DSMC  computations  are  begun  at  a  location  downstream  of 
the  nozzle  throat.  The  initial  profiles  of  flow  properties  are  obtained  from  continuum  simulations  provided 
by  the  code  of  Butler.^  The  interaction  of  the  gas  with  the  wall  of  the  nozzle  is  assumed  to  be  diffuse  with  full 
energy  accommodation  to  the  wall  temperature  profile  obtained  from  the  continuum  simulation.  Charged 
species  are  neutralized  upon  contact  with  the  wall.  The  computations  extend  through  the  nozzle  exit  plane 
to  a  distance  of  30  mm  from  the  thruster  along  the  axis.  Expansion  is  assumed  to  occur  into  a  perfect 
vacuum. 

The  grids  employed  typically  consist  of  500  by  55  cells.  Due  to  the  wide  dynamic  range  in  density  for 
these  flows,  the  cells  are  adapted  to  the  local  mean  free  path,  and  the  time  step  in  each  cell  is  scaled  by  the 
local  mean  time  between  collisions.  Despite  the  time  step  variation,  the  transient  phase  of  the  simulation 
remains  long.  Typically,  it  requires  at  least  15,000  iterations  before  steady  state  is  reached.  Macroscopic 
results  are  then  obtained  by  sampling  properties  over  a  further  10,000  time  steps.  At  steady  state,  there  are 
at  least  300,000  particles  present  in  the  simulation. 

Results 

Results  are  presented  for  a  mass  flow  rate  of  13.3  mg/sec  and  power  levels  of  1.34  and  0.80  kW.  These 
correspond  very  closely  to  conditions  considered  in  a  number  of  experimental  investigations.  Comparison  of 
DSMC  predictions  with  experimental  data  for  various  flow  field  properties  is  made  at  the  following  locations: 
(1)  in  the  interior  of  the  nozzle;  (2)  in  the  nozzle  exit  plane;  and  (3)  in  the  plume  flow  field.  In  addition, 
comparison  is  made  between  measurement  and  prediction  for  specific  impulse. 
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Nozzle  Interior 


In  Fig.  2a,  the  variation  of  axial  velocity  component  along  the  nozzle  centerline  is  shown  at  a  power  of 
1.34  kW.  The  axis  labeling  is  such  that  an  axial  coordinate  of  zero  is  located  at  the  nozzle  exit  plane.  The 
experimental  data  is  due  to  Storm  and  Cappelli.^^  Note  that  these  data  were  measured  at  a  slightly  higher 
arcjet  power  of  1.48  kW.  The  DSMC  simulation  provides  quite  good  agreement  with  the  data  although  there 
is  a  tendency  to  predict  slightly  lower  values.  This  may  partly  be  attributed  to  the  lower  power  considered  in 
the  simulation.  Far  upstream,  the  data  indicates  velocities  that  are  significantly  higher  than  the  computed 
values;  however,  this  is  the  region  of  greatest  experimental  uncertainty.  The  minimum  turning  point  shown 
in  the  DSMC  calculations  is  a  result  of  the  diffusion  of  cold  hydrogen  molecules  from  the  boundary  layer 
into  the  hot  core  region  where  the  plasma  deposition  of  energy  is  greatest. 

Radial  profiles  of  axial  velocity  component  at  three  different  axial  locations  inside  the  nozzle  are  shown 
in  Fig.  2b.  In  general,  the  shapes  of  the  measured  and  predicted  profiles  are  in  good  agreement  and  are 
found  to  improve  with  increased  downstream  distance.  In  particular,  the  profile  computed  at  a  distance  of 
just  -2.54  mm  from  the  nozzle  exit  plane  shows  the  best  agreement  with  the  data. 

Nozzle  Exit  Plane 

In  Fig.  3a,  the  concentrations  of  molecular  and  atomic  hydrogen  across  the  nozzle  exit  plane  are  shown. 
Experimental  measurement  of  H2  was  performed  by  Beattie  and  Cappelli^^  using  Raman  scattering.  The 
measurements  for  H  are  those  reported  by  Wysong  et  al.^^  obtained  using  a  two-photon  laser  induced 
fluorescence  technique.  In  general,  the  agreement  between  computation  and  measurement  is  better  for  H2 
than  for  H.  In  particular,  the  simulations  predict  an  atomic  concentration  on  the  axis  that  is  a  factor  of 
3  to  4  smaller  than  the  experimental  data.  A  further  simulation  performed  with  the  dissociation  rate  for 
H2“H  collisions  increased  by  a  factor  of  10  led  to  an  increase  in  H-atom  concentration  at  the  nozzle  exit 
of  just  10%.  The  discrepancy  between  measurement  and  calculation  for  the  hydrogen  atom  concentrations 
is  therefore  attributed  either  to  too  little  dissociation  in  the  continuum  solution  used  to  start  the  DSMC 
calculation,  or  to  experimental  error. 

In  Fig.  3b,  the  translational  temperature  profile  of  atomic  hydrogen  across  the  nozzle  exit  plane  is 
shown.  Two  sets  of  experimental  data  reported  by  Wysong  et  al.^^  and  by  Storm  and  Cappelli^^  are  shown. 
In  general,  the  DSMC  results  show  a  small  underprediction  of  the  temperature,  although  the  agreement  is 
quite  good.  The  degree  of  thermal  nonequilibrium  of  the  flow  in  the  nozzle  exit  is  illustrated  in  Fig.  3c. 
Here,  the  translational,  rotational,  and  vibrational  temperatures  of  molecular  hydrogen  predicted  by  the 
DSMC  technique  are  compared  with  the  rotational  temperature  measurements  of  Beattie  and  Cappelli.^^ 
The  predicted  rotational  temperatures  offer  satisfactory  agreement  with  the  measured  data.  The  peak 
vibrational  temperature  is  about  1,000  K  higher  than  the  peak  translational  temperature. 

A  further  indication  of  the  degree  of  thermal  nonequilibrium  in  this  flow  is  illustrated  in  Fig.  3d  in 
which  the  electron  temperature  at  the  nozzle  exit  plane  is  shown.  The  experimental  data  is  due  to  Storm 
and  Cappelli.^^  There  are  two  DSMC  results  shown  in  which  the  effect  of  including  ohmic  dissipation 
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(r^)  is  considered.  With  ohmic  dissipation,  there  is  quite  good  agreement  obtained  with  the  experimental 
data.  By  comparison,  omission  of  H  leads  to  significant  underprediction  of  the  electron  temperature.  Note 
that  with  ohmic  heating,  the  peak  electron  temperature  is  an  order  of  magnitude  higher  than  the  peak 
translational  temperature,  again  illustrating  the  strong  degree  of  thermal  nonequilibrium  in  these  flows.  In 
Fig.  3e,  the  electron  concentration  across  the  nozzle  exit  plane  is  shown.  The  experimental  data  was  also 
reported  in  Ref.  24.  The  agreement  between  simulation  and  measurement  is  very  poor  with  the  DSMC 
technique  predicting  concentrations  that  are  an  order  of  magnitude  higher  than  the  data.  This  behavior  is 
attributed  to  the  very  simple  model  employed  in  the  DSMC  computation  for  simulating  the  electric  field 
effects.  Remember,  the  electrons  are  constrained  to  move  with  the  average  ion  velocity,  and  this  is  obviously 
much  slower  than  would  occur  in  the  actual  flow. 

Finally,  the  axial  component  of  velocity  in  the  nozzle  exit  plane  is  shown  in  Fig.  3f.  The  experimental 
data  is  that  of  Wysong  et  al..^^  While  offering  relatively  good  agreement,  the  DSMC  results  are  consistently 
lower  than  the  measured  data  across  most  of  the  exit  plane.  The  shape  of  the  velocity  profile  is  captured 
very  well  by  the  simulation. 

In  Figs.  4a-4d,  comparisons  of  simulation  and  data  for  a  power  level  of  0.80  kW  are  shown  for  neutral 
species  concentrations,  translational  temperature  of  hydrogen  atoms,  internal  temperatures  of  hydrogen 
molecules,  and  axial  velocity  component  of  hydrogen  atoms.  Overall,  the  same  observations  hold  at  this 
power  level  as  for  the  1.34  kW  case.  The  agreement  between  simulation  and  measurement  is  generally  very 
good  except  for  the  hydrogen  atom  concentrations  in  which  the  DSMC  technique  predicts  densities  that  are 
a  factor  of  about  3  lower  than  the  experimental  data. 

Plume  Flow  Field 

Of  the  many  experimental  investigations  of  hydrogen  arcjet  flows.  Ref.  22  is  unique  in  that  data  taken 
in  the  plume  are  reported  that  are  unaffected  by  chamber  effects  are  reported.  These  data  are  compared  in 
this  section  with  the  corresponding  DSMC  results  out  to  a  distance  of  30  mm  from  the  nozzle  exit  plane. 
Variation  of  translational  temperature  and  axial  velocity  component  for  atomic  hydrogen  along  the  axis  of 
the  plume  are  shown  in  Figs.  5a  and  5b  for  a  power  level  of  1.34  kW.  In  these  plots,  the  axis  labeling  is  such 
that  the  nozzle  exit  plane  is  again  located  at  zero.  It  should  be  noted  that  these  are  the  first  comparisons 
between  simulation  and  measurement  conducted  in  the  plume  of  an  arcjet  thruster.  As  shown  in  Fig.  5a, 
the  temperature  profile  predicted  by  the  DSMC  technique  is  in  excellent  agreement  with  the  experimental 
data.  For  velocity,  the  DSMC  results  lie  consistently  1-2  km/s  below  the  data;  however,  it  is  significant  that 
the  lack  of  change  in  velocity  is  predicted  by  the  simulation. 

Axial  plume  profiles  of  translational  temperature  and  axial  velocity  component  at  a  power  level  of 
0.80  kW  are  shown  in  Figs.  6a  and  6b.  In  this  case,  the  DSMC  temperature  predictions  are  distinctly  lower 
than  the  measured  data.  The  velocity  predictions  are  again  lower  than  the  experimental  results  by  1-2  km/s. 

In  the  experimental  investigation  of  Wysong  et  al.^^,  radial  profiles  of  plume  properties  were  also  mea¬ 
sured  at  an  axial  distance  of  10  mm  from  the  nozzle  exit  plane.  In  Figs.  7a  and  7b  profiles  of  translational 
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temperature  and  axial  velocity  component  of  atomic  hydrogen  are  shown  for  a  power  of  1.34  kW.  Within 
the  large  scatter  exhibited  by  the  measurements,  the  peaked  profile  predicted  by  the  DSMC  technique  offers 
reasonable  agreement.  The  velocity  profiles  are  in  excellent  agreement. 

In  Figs.  8a  and  8b,  the  temperature  and  velocity  profiles  for  a  power  of  0.80  kW  are  shown.  The  DSMC 
method  predicts  significantly  lower  temperatures  and  slightly  lower  velocities. 

Axial  and  radial  profiles  of  hydrogen  atom  concentration  were  also  measured  in  Ref.  22.  The  comparisons 
of  those  data  with  the  DSMC  predictions  are  consistent  with  the  behavior  at  the  nozzle  exit  plane.  Thus, 
in  a  fairly  consistent  way,  the  DSMC  method  predicts  densities  for  hydrogen  atoms  throughout  the  plume 
that  are  a  factor  of  3  to  4  lower  than  the  measured  data. 

Thruster  Performance 

Validation  of  DSMC  predictions  for  fiow  field  properties  is  important  for  assessment  of  spacecraft  in¬ 
teraction  effects  such  as  contamination.  However,  an  equally  important  role  for  numerical  simulations  lies 
in  the  optimization  of  arcjet  design.  For  computational  methods  to  be  taken  seriously  in  this  respect,  it  is 
essential  that  they  are  able  to  accurately  predict  overall  thruster  performance.  A  series  of  measurements 
of  the  thrust  obtained  from  the  hydrogen  arcjet  under  different  operating  conditions  and  geometries  was 
conducted  several  years  ago  at  NASA  Lewis  Research  Center.^  The  data  employed  here  for  comparison  with 
the  DSMC  results  is  for  the  arcjet  configuration  that  most  closely  resembled  that  employed  in  the  flow  field 
measurements  and  simulations.  In  Fig.  9,  variation  of  specific  impulse  with  specific  energy  (power/fiow  rate) 
is  shown.  Experimental  data  is  provided  at  two  different  flow  rates  that  bound  the  condition  considered  in 
the  simulations.  Two  sets  of  DSMC  data  are  shown  that  indicate  the  effect  of  including  ohmic  dissipation  in 
the  simulation.  As  expected,  with  finite  H,  the  specific  impulse  is  increased,  and  the  amount  of  the  increase 
appears  to  grow  with  specific  energy.  Due  to  the  scatter  in  the  measurements  and  the  fact  that  the  flow  rates 
are  different  from  the  simulation,  it  is  difficult  to  make  a  definitive  comparison.  However,  it  does  appear 
that  the  DSMC  predictions  that  include  ohmic  heating  agree  with  the  measured  data  to  within  5%.  Further 
computations  need  to  be  undertaken  to  cover  a  wider  range  of  operating  conditions. 

Concluding  Remarks 

A  Monte  Carlo  method  has  been  developed  and  applied  to  simulate  flow  in  the  nozzle  and  plume  of  a 
hydrogen  arcjet.  The  model  includes  detailed  simulation  of  nonequilibrium  collision  and  plasma  phenomena. 
These  models  are  implemented  in  a  numerically  efficient  computer  code.  In  this  study,  extensive  comparison 
between  code  predictions  and  a  variety  of  experimental  measurements  were  presented.  The  comparisons 
were  made  for  two  different  arcjet  power  levels,  and  for  several  different  aspects  of  the  arcjet  flow  fields. 
In  particular,  predicted  results  were  compared  with  measured  data  in  the  nozzle  interior,  at  the  nozzle  exit 
plane,  and  in  the  plume  expansion.  Monte  Carlo  prediction  of  specific  impulse  was  also  compared  with 
experimental  measurement. 

In  terms  of  the  arcjet  flow  fields,  there  was  generally  good  agreement  between  prediction  and  measure- 
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ment  for  both  translational  temperature  and  velocity.  There  was  an  overall  tendency  in  the  simulations  to 
predict  slightly  lower  values  of  temperature  and  velocity  in  comparison  with  the  experimental  results.  The 
greatest  discrepancy  occurred  for  the  concentration  of  atomic  hydrogen  for  which  the  Monte  Carlo  approach 
predicted  values  that  were  a  factor  of  3  to  4  lower  than  those  measured.  This  behavior  was  found  to  be 
consistent  at  the  nozzle  exit  and  throughout  the  plume  flow  field.  A  significant  increase  in  the  hydrogen 
dissociation  rate  employed  in  the  simulation  did  not  increase  the  hydrogen  atom  concentration  significantly. 
Therefore,  it  is  postulated  that  the  discrepancy  lies  either  in  the  starting  conditions  employed  in  the  Monte 
Carlo  calculations,  or  in  the  calibration  of  the  diagnostic  technique. 

At  the  nozzle  exit  plane,  the  nozzle  expansion  process  has  created  a  flow  that  is  in  a  state  of  strong 
thermal  nonequilibrium.  The  Monte  Carlo  simulation  predicts  different  values  of  temperature  for  the  trans¬ 
lational,  rotational,  and  vibrational  energy  modes,  and  for  the  electrons.  In  terms  of  the  electrons,  it  is 
demonstrated  that  a  model  is  required  to  simulate  ohmic  heating  in  order  to  accurately  predict  the  electron 
temperature  profiles  measured  experimentally.  It  was  also  found  that  a  simple  model  employed  to  simulate 
electric  field  effects  results  in  prediction  of  electron  concentrations  that  were  too  high  by  about  an  order  of 
magnitude. 

For  thruster  performance,  the  Monte  Carlo  simulation  is  able  to  reproduce  experimental  measurements 
within  about  5%  when  ohmic  heating  is  included  in  the  computation. 

In  summary,  the  generally  successful  comparisons  between  prediction  and  measurement  of  flow  field 
properties  indicates  that  the  Monte  Carlo  method  is  capable  of  computing  hydrogen  arcjet  flows  with  en¬ 
gineering  accuracy.  This  has  the  significance  of  demonstrating  that  the  numerical  method  may  be  used  for 
prediction  of  spacecraft  interaction  effects  when  arcjets  are  employed  on  real  satellite  configurations.  The 
good  agreement  obtained  between  computation  and  measurement  for  specific  impulse  demonstrates  that  the 
Monte  Carlo  method  may  be  employed  for  optimization  of  the  performance  of  hydrogen  arcjets. 

Additional  research  on  the  treatment  in  the  Monte  Carlo  approach  of  electric  field  effects  should  lead  to 
even  better  correspondence  between  prediction  and  measurement  for  both  flow  field  properties  and  thruster 
performance.  For  application  of  the  Monte  Carlo  method  to  the  hydrazine  arcjets  that  are  flown  in  space, 
work  is  also  needed  to  develop  appropriate  models  for  this  more  complicated  thermochemical  nonequilibrium 
system. 
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Fig.  1.  Schematic  diagram  of  the  1  kW  hydrogen  arcjet. 
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Fig.  2a.  Axial  velocity  component  along  the  nozzle  centerline  for  P=1.34  kW. 


Fig.  2b.  Radial  profiles  of  axial  velocity  component  at  different  axial  locations  in  the  nozzle  for  P=1.34  kW. 
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Fig.  3b.  Translational  temperature  of  atomic  hydrogen  in  the  nozzle  exit  plane  for  P=1.34  kW. 
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Fig.  3c.  Thermal  nonequilibrium  of  molecular  hydrogen  in  the  nozzle  exit  plane  for  P=1.34  kW 


Fig.  3d.  Electron  temperature  in  the  nozzle  exit  plane  for  P=1.34  kW 
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